Passive incoherent imaging techniques such as phase diversity and wavelength diversity have long been successfully used to compensate for atmospheric turbulence. Both diversity schemes traditionally use iterative 2-D Fourier transforms in a sequential error-reduction methodology that is typically slow and relegated to post-processing applications.
For optical imaging systems having apertures that are significantly larger than the atmospheric coherence length, ro, the turbulent atmosphere usually is the most significant contributor to the imaging system's loss of resolution [1]. Atmospheric turbulence introduces aberrations that can significantly degrade the performance of an optical imaging system. The degradation in imaging performance depends on a variety of factors, such as the operating wavelength(s), the relative size of the imaging system's entrance pupil diameter as compared to the atmospheric coherence length, ro (Fried parameter), the angular instantaneous field of view (IFOV) as compared to the isoplanatic angle, θ0, sampling effects, signal-to-noise issues, and system effects such as imperfect optics, fixed pattern noise and platform jitter.
Given a circular aperture, the upper bound on achievable resolution is the so-called diffraction limit,
                                          Δ            ⁢                                                  ⁢                          x              diff                                =                      1.22            ⁢                                                  ⁢                                          λ                _                            D                        ⁢            z                          ,                            (        1        )            where λ is the center wavelength of the illuminating light, D is the diameter of the entrance pupil of the imaging system (for example, the diameter of the telescope, or camera lens), and z is the distance between the imaging system's entrance pupil and the object to be imaged. For example, an object located 2 km from a camera with a 4½ inch lens would have a diffraction-limited resolution of 1.059 cm at a center wavelength of 500 nm. This means that object features less than this size will not be resolved and object features larger than this size will be resolved if diffraction-limiting imaging conditions are attained.
Unfortunately, the atmosphere severely degrades the ability to achieve the diffraction-limited resolution indicated in equation (1). Instead, the conventionally-attainable maximum resolution of an imaging system that is looking through atmospheric turbulence is given by
                                          Δ            ⁢                                                  ⁢                          x                              at                ⁢                                                                  ⁢                m                                              =                      1.22            ⁢                                                  ⁢                                          λ                _                                            r                0                                      ⁢            z                          ,                            (        2        )            where the expression in the denominator of equation (2) is the atmospheric coherence length, or Fried parameter. For imaging over horizontal paths, recent experimental values for ro have ranged from 1 to 4 cm leading to a “loss” of resolution of 11.52 (for an ro of 1 cm) or 2.88 (for an ro of 4 cm) with respect to the diffraction-limited resolution given in equation (I). Stated another way, by compensating for the effects of atmospheric turbulence, a maximum increase in resolution between 2.88 and 11.52 can be expected for the center wavelength and entrance-pupil diameter specified above with a suitable turbulence-compensation approach. These effects are even more dramatic when a telescopic imaging system is used. For instance, if an eight-inch telescope is attached to the camera or video-camera at optical wavelengths, the maximum increase in spatial resolution for the same values of ro above jumps to 20.32 (for an ro of 1 cm) or 5.08 (for an ro of 4 cm).
As can be seen by dividing equation (2) by equation (1), that the expected increase in resolution in a linear direction of an imaging system with entrance-pupil diameter D is given by
                    R        =                                            Δ              ⁢                                                          ⁢                              x                                  a                  ⁢                                                                          ⁢                  tm                                                                    Δ              ⁢                                                          ⁢                              x                diff                                              =                                    D                              r                0                                      .                                              (        3        )            
With knowledge of the diameter of the imaging system's entrance pupil and the value of ro for an illuminating wavelength, equation (3) can be used to determine the maximum achievable increase in resolution, neglecting system effects and assuming full compensation of atmospheric turbulence.
If the value of ro is known at a given wavelength, then it can be scaled to another wavelength by
                                          r            0                          λ              2                                =                                    r              0                              λ                1                                      ·                                          (                                                      λ                    2                                                        λ                    1                                                  )                                            6                5                                                    ,                            (        4        )            where r0λ1 is the value of ro at λ1 and r0λ2 is the value of ro at wavelength of interest λ2. For example, if the value of ro at the illuminating wavelength of 500 nm is 1 cm, then the value of ro at 1.06 μm is about 2.46 cm. Once the new value of ro is known, then equation (3) can be used to determine the upper bound on the resolution increase of the optical system under consideration.
Equation (4) can be obtained from
                                          r            0                    =                      0.185            ⁢                                          (                                                      λ                    2                                                                              ∫                      0                      z                                        ⁢                                                                                            c                          n                          2                                                ⁡                                                  (                          ξ                          )                                                                    ⁢                                              ⅆ                        ξ                                                                                            )                                            3                5                                                    ,                            (        5        )            where the term in the denominator is the integral of the atmospheric structure constant along the optical path [0, z].
In order to attain an increase in resolution approaching the limits given by equation (3), the effects of atmospheric turbulence must be deconvolved from the aberrated image. FIG. 1 shows the effects of simulated atmospheric turbulence for a target that is 1.3 km away from the optical imaging system. The center wavelength was 550 nm and the distance between the simulated object and the 8-inch telescopic imaging system was 1.3 km. The Fried coherence length was 2 cm. The imaging system itself was assumed to be noise-free to isolate and demonstrate the atmospheric turbulence effects on the resulting image.
FIG. 2 shows the simulated diffraction-limited, atmospheric-turbulence-compensated image from the same imaging system. Notice the dramatic increase in spatial resolution possible by compensating for the effects of atmospheric turbulence. From FIG. 1, it can be seen that even if a “perfect” (aberration-free) imaging system could be built (for example, no optical system aberrations), the presence of atmospheric turbulence still severely degrades the imaging system's performance in terms of achievable spatial resolution and that compensating or correcting for the effects of atmospheric turbulence can significantly improve the overall image quality.
For comparisons to real data results, FIG. 3 shows actual imagery taken under imaging conditions similar to those that were used to generate the simulated images above. The data in FIG. 3 was taken by Dr. C. J. Carrano from Lawrence Livermore National Labs (LLNL) using a very promising speckle imaging-based methodology [2]. Dr. Carrano implements her turbulence-compensation method on a Field Programmable Gated Array (FPGA) technology and the image-processing time of her method is expected to perform near real-time [3,4]. Note the similarity in image quality between the simulated and real-image data. The process of removing atmospheric turbulence from images is called deconvolution.
There are three basic ways to estimate and remove the effects of atmospheric turbulence, namely, 1) adaptive optics systems, 2) post-processing atmospheric turbulence compensation systems, and 3) hybrid systems.
Adaptive optics systems are hardware-based systems that can correct atmospheric-turbulence effects in real-time (faster than 30 Hz). Adaptive optics systems are generally cumbersome, require extensive hardware, alignment, and expertise, and are expensive. They are also predominantly designed for fixed sites and aren't typically man-portable, rugged, or covert.
Post-processing atmospheric turbulence compensation systems are largely implemented in software (software dominant) but traditionally are very slow (not real-time). Hybrid methods are a cross between the two and generally relax some of the hardware processing requirements and then correct for the loss in performance using post-processing methods.
Many of the existing software-dominant deconvolution approaches such as wavelength diversity, phase diversity, multi-frame blind deconvolution and other post-processing approaches have been predominantly focused on phase-only corrections (near-field turbulence approximation) and have been iterative and slow[5,6,7]. Some notable exceptions are Paxman's use of phase diversity for reconstructing the unaberrated object brightness through distributed turbulence, Aubailly and Vorontsov's methods for fusing segments of “lucky” image frames to de-blur images obtained over horizontal optical paths, and more recently the work of Carrano which generalized the traditional speckle-imaging technique that is and has been used to great effect in astronomy and other near-field turbulence application areas [3,4,8,9]. The approach of Ortiz and Carrano et al. improves the speed of the traditional speckle imaging-based post-processing Methods. For this reason a quick overview of Carrano's method is provided along with some relevant comparisons and observations.
Carrano's approach parameterized the Korff transfer function in terms of the atmospheric coherence length, ro to estimate the object magnitude and used the bi-spectrum technique to estimate the object phase without requiring the presence of a reference source such as a star or laser guide-star [10]. The approach was applied to a horizontal path turbulence scenario with excellent results for optical path lengths of 0.5 km to 10 km. Processing was done on a 256 by 256 pixel by 100 image frame data cube using a 1.7 GHz Pentium IV processor. It originally took about 10 seconds to process this data cube, which represented an isoplanatic patch from a 1280 by 1024 image. To recover the full image, 20 similar data cubes must be processed resulting in 200 seconds of processing time for their non-optimized algorithm. In subsequent papers, Ortiz, Carrano et al. have continued to improve processing times to approximately 1 second using faster processing capabilities [3,4]. When considering real-time atmospheric turbulence compensating systems that work with uncooperative targets (for example, in some surveillance applications), a distinct limitation of the Carrano approach is that it could require hundreds of images or more to provide the parameters needed for their method to work. It would be advantageous to require just one pair of simultaneously-captured images given sufficient signal-to-noise to form the image pair (for example, enough signal to capture the blurry pair of images).
The following is an overview of the optical-systems model used in general-purpose incoherent optical imaging systems. This will provide a common basis of understanding when comparing alternative imaging methodologies and will provide the necessary background to understand the approach of the present invention.
For many incoherent imaging applications, a linear, shift-invariant imaging model is appropriate. For this case, the optical systems model is given by [11],i({right arrow over (x)})=o({right arrow over (x)})*|hi({right arrow over (x)})|2,  (6)where o({right arrow over (x)}) the 2-D pristine object brightness function, |hi({right arrow over (x)})|2 is the imaging system's point spread function (PSF), i({right arrow over (x)}) is the “blurry” image due to atmosphere and optical imaging system effects, and {right arrow over (x)} is a 2-D position vector in the image plane. The asterisk represents 2-D spatial convolution.
FIG. 4 shows the general relationship among the quantities of interest in the incoherent imaging scenario in a general-purpose optical systems model. By taking the 2-D Fourier transform of both sides of equation (6), the frequency-space equivalent of equation (6) is given byI({right arrow over (f)})=O({right arrow over (f)})H({right arrow over (f)}),  (7)where I({right arrow over (f)}) is the image spectrum, O({right arrow over (f)}) is the object spectrum, H({right arrow over (f)}) is the optical transfer function (OTF) and {right arrow over (f)} is a 2-D spatial frequency variable. Equations (6) and (7) apply at each spatial coordinate {right arrow over (x)} and at each spatial frequency {right arrow over (f)}.
The PSF can be related to the optical transfer function by,
                                          H            ⁡                          (                              f                →                            )                                =                                                  ⁢                              ⌊                                                                                                                        h                        i                                            ⁡                                              (                                                  x                          →                                                )                                                                                                  2                                ⌋                                                                                    ⁡                                  [                                                                                                                                    h                          i                                                ⁡                                                  (                                                      x                            →                                                    )                                                                                                            2                                    ]                                                                              f                  _                                =                0                                                    ,                            (        8        )            where the symbol, ℑ[•], denotes taking the 2-D Fourier transform of the expression inside the brackets. The optical transfer function is seen to be the 2-D Fourier transform of the PSF and then normalized so that the maximum value of the OTF is 1.
The relationship between the optical system's impulse response hi({right arrow over (x)}) and the generalized pupil function (GPF) is given by,hi({right arrow over (x)})=ℑ−1[W({right arrow over (x)})],  (9)where ℑ−1[•] represents taking the 2-D inverse Fourier transform of the expression inside the brackets, andW({right arrow over (x)})=A({right arrow over (x)})ejφ({right arrow over (x)}),  (10)is the generalized pupil function (GPF). The function A({right arrow over (x)}) is an aperture function that has a value of 1 inside the clear aperture of the telescope and 0 outside of it. The function φ({right arrow over (x)}) is the atmospherically-induced phase aberration realized at spatial coordinate {right arrow over (x)} in the entrance pupil of the imaging system. In the near-field turbulence approximation model, such as if one is looking up through the atmosphere into space from a ground-based telescope, the amplitude variations are assumed negligible within an isoplanatic patch, and are set to 1.
In traditional “diversity-based” post-processing atmospheric turbulence compensation methods, the basic idea is to insert a known diversity into the imaging system by some artifice and then simultaneously capture the original image and the diversity image. The OTFs of the aberrated image and the diversity image are related to each other and a suitable error metric is used to select the OTF that produces the lowest local and global error. Depending on the approach taken, in relating the OTF to the diversity OTF, often the resulting error metric will be a function of the aberrated image spectrum, the diversity image spectrum (both measured) and analytical expressions of the OTF and diversity OTF that are both functions of entrance pupil phase estimates or phase difference estimates. For instance, the popular phase-diversity post-processing atmospheric turbulence compensation method uses an additive phase term that is known a priori in the expression of the phase diversity generalized pupil function,Wpd({right arrow over (x)})=A({right arrow over (x)})ej(φ({right arrow over (x)})+φpd({right arrow over (x)})),  (11)where the subscript pd denotes the phase diversity method was used. The expression Wpd({right arrow over (x)}) is the phase diversity generalized pupil function and, as shown, has a known phase diversity φpd({right arrow over (x)}) added to the unknown atmospheric turbulence phase φ({right arrow over (x)}) at every entrance pupil spatial coordinate {right arrow over (x)}. Often a known quadratic phase factor can be introduced in the phase diversity image by slightly defocusing the diversity image.
Another diversity method is one by which an image is captured simultaneously at two different narrow-band wavelengths centered at λ1 and λ2. The wavelength diversity generalized pupil function is then given by
                                          W            wd                    ⁡                      (                          x              →                        )                          =                              A            ⁡                          (                              x                →                            )                                ⁢                                    ⅇ                              j                ⁡                                  (                                                                                    λ                        1                                                                    λ                        2                                                              ⁢                                          ϕ                      ⁡                                              (                                                  x                          →                                                )                                                                              )                                                      .                                              (        12        )            In the traditional diversity-based atmospheric turbulence compensation methods, a diversity OTF is generated by using equations (8) through (10) in reverse order and substituting the appropriate diversity generalized pupil function from equation (11) or (12) depending on which diversity method one is using for equation (10). A common error metric such as the Gonsalvez error metric [12,13],
                              E          ⁡                      (                          f              →                        )                          =                                                                                                                I                    ⁡                                          (                                              f                        →                                            )                                                        ⁢                                                                                    H                        ^                                            d                                        ⁡                                          (                                              f                        →                                            )                                                                      -                                                                            I                      d                                        ⁡                                          (                                              f                        →                                            )                                                        ⁢                                                            H                      ^                                        ⁡                                          (                                              f                        →                                            )                                                                                                          2                                                                                                                    H                    ^                                    ⁡                                      (                                          f                      →                                        )                                                                              2                        +                                                                                                                      H                      ^                                        d                                    ⁡                                      (                                          f                      →                                        )                                                                              2                                                          (        13        )            can then be applied at every point in the image spectrum as a means to determine when the OTF estimate is accurate enough. Note that the carat symbol ^ above the OTF and diversity OTF indicate that these quantities are estimated. In traditional diversity-based atmospheric turbulence compensation methods, the process for estimating the OTF (and also the diversity OTF by analogy) includes                1. Using a suitable basis set like the Zernike polynomials initially to generate an entrance pupil plane-phase estimate. This is done by just initially guessing the phase—for example, that all phase values are zero;        2. Forming the generalized pupil function with this entrance pupil phase “guess” using equation (10). Note: use equations (11) or (12) for estimating the diversity OTF;        3. Zero-packing the GPF for sampling reasons in preparation of generating an OTF estimate;        4. Forming the impulse response in accordance with equation (9);        5. Determining the PSF estimate by point-wise taking the magnitude squared of the result of step 4 above;        6. Using equation (8) to form the OTF estimate;        7. After forming both the OTF and diversity OTF estimate, applying an error metric such as the Gonsalvez error metric given in equation (13) and possibly some constraints to determine the instantaneous error at each spatial location in the OTF;        8. Summing the errors to determine the total summed error due to the initial entrance pupil plane phase estimation error;        9. Changing the weights on the Zernike polynomials in a methodical manner to come up with a new entrance pupil plane phase estimate;        10. Repeating steps 2 through 8 to generate a new error estimate;        11. Comparing the new error estimate to the old estimate and keeping the phase estimates associated with the lowest error; and        12. Continuing to execute steps 10 and 11 until the error is minimized and the best OTF estimate is obtained.        
Once the error-minimized OTF estimate is obtained, a Wiener filter can be generated that removes the effect of atmospheric turbulence. Notice that the OTF itself is generated by phase estimates that are due to 1) atmospheric turbulence effects, and 2) aperture effects (for example, diffraction effects). If the effects of the atmospheric turbulence are mitigated by filtering them out using the Wiener filter, then the only remaining effect is that due to diffraction, and so the diffraction-limited result is obtained. To attempt to remove the effects of the aperture, super-resolution methods need to be employed.
After the error-minimized OTF estimate is determined, the Wiener filter is given by
                                                        H                              -                1                                      ⁡                          (                              f                →                            )                                =                                                    H                *                            ⁡                              (                                  f                  →                                )                                                    (                                                                                                              H                      ⁡                                              (                                                  f                          →                                                )                                                                                                  2                                +                α                            )                                      ,                            (        14        )            where the asterisk on the right side of equation (14) represents complex conjugation. Care must be taken for the case where the denominator of equation (14) approaches zero. A parameter α based on system noise is sometimes included in the denominator to prevent equation (14) from blowing up as H({right arrow over (f)}) approaches zero. As can be seen from equation (7), multiplying the image spectrum by the Wiener filter leads to an unaberrated object spectrum,O({right arrow over (f)})=I({right arrow over (f)})H−1({right arrow over (f)}),  (15)and the atmospheric turbulence free object brightness estimate is simply obtained by taking the 2-D inverse Fourier transform of equation (15).
Notice also that equation (8) can be directly determined from the GPF in the following manner,
                              H          ⁡                      (                          f              →                        )                          =                                                            W                ⁡                                  (                                      x                    →                                    )                                            ⊗                              W                ⁡                                  (                                      x                    →                                    )                                                                                    W                ⁡                                  (                  0                  )                                            ⊗                              W                ⁡                                  (                  0                  )                                                              .                                    (        16        )            The symbol  means auto-correlation and the entrance pupil spatial position variable is related to the spatial frequency variable by{right arrow over (x)}=λdi{right arrow over (f)},  (17)where di is the distance from the imaging system's exit pupil to the focal plane and lambda is the center wavelength of the illuminating light. The denominator in equation (16) is just the area of the imaging system's entrance pupil (commonly the area of the telescope's collecting lens for well-designed optical systems).
In the past, due to limitations of processor speed and some computational inefficiencies inherent in the majority of traditional post-processing atmospheric turbulence compensation approaches, this direct methodology illustrated by equation (16) was considered from a practical point of view to be much slower than the iterative Fourier transform-based approach described in steps 1 through 12 above. What is needed is a new methodology that overcomes these inefficiencies and provides for a software-dominant approach for atmospheric turbulence compensation that can be accomplished in real-time. It would be beneficial to both significantly modify the traditional atmospheric turbulence compensation methodology and also use a general-purpose parallel-processing device such as a field-programmable gated array to achieve real-time processing speeds. Utilizing both steps would enable a system to achieve real-time processing in a practical and scalable manner.